package edu.vt.cbil.core;

import edu.vt.cbil.dataset.*;
import edu.vt.cbil.util.BioToolbox;
import edu.vt.cbil.util.ChartXYLine;
import edu.vt.cbil.util.ConfigObject;
import edu.vt.cbil.util.Histogram;
import org.apache.commons.math3.stat.correlation.PearsonsCorrelation;
import org.apache.commons.math3.transform.DctNormalization;
import org.apache.commons.math3.transform.FastCosineTransformer;
import org.apache.commons.math3.transform.TransformType;
import org.jfree.chart.plot.XYPlot;
import org.jfree.chart.renderer.xy.XYLineAndShapeRenderer;
import org.jfree.data.xy.XYSeries;

import java.awt.*;
import java.io.File;
import java.text.SimpleDateFormat;
import java.util.*;
import java.util.List;
import java.util.logging.Logger;

/**
 * The type Segment processor.
 *
 * @author Justin
 * @version 1.0
 * @created 12 -11-2014 21:40:20
 */
public class SegmentProcessor {
    /**
     * The M _ cutter.
     */
    public CNSegmentCutter m_cutter;
    private Logger m_log;
    private BioCollection p_Collection;
    private ConfigObject p_config;
    private String curSampleId;
    private BioSample p_sample;
    private double minCN_segment;
    private TreeMap<Integer, Double> CNPeakArray;
    private TreeMap<Double, Boolean> LOHArray;

    /**
     * Instantiates a new Segment processor.
     *
     * @param pc the pc
     */
    public SegmentProcessor(BioCollection pc, ConfigObject pcfg) {
        m_log = Logger.getLogger("edu.vt.cbil.bacom2");

        p_Collection = pc;
        p_config = pcfg;
        m_cutter = new CNSegmentCutter(pc,
                p_config.params.PVALUE_THRESH, p_config.params.MIN_SEG_LENGTH,
                p_config.params.MIN_SEG_STEP, p_config.params.COMBINE_FACTOR);
        CNPeakArray = null;
        LOHArray = null;
    }

    /**
     * Sets collection.
     *
     * @param pc the pc
     */
    public void setCollection(BioCollection pc) {
        this.p_Collection = pc;
        m_cutter.setCollection(pc);
    }

    /**
     * Sets current sample.
     *
     * @param sampleId the sample id
     */
    public void setCurrentSample(String sampleId) {
        curSampleId = sampleId;
        p_sample = p_Collection.getSample(sampleId);
        m_cutter.setSample(p_sample);

    }

    private void smoothingCopyNumber(int width) {
        m_cutter.prepareCopyNumberSegment(p_sample);
        if (width > 2) {
            m_log.fine("Smoothing Width = " + width);
            m_cutter.smoothing(width);
        }

    }

    /**
     * Cut segment.
     *
     * @param sampleId    the sample id
     * @param smoothWidth the smooth width
     */
    public void cutSegment(String sampleId, int smoothWidth) {
        setCurrentSample(sampleId);

        smoothingCopyNumber(smoothWidth);

        m_cutter.splitChromosomes(p_sample);

        clearCutter();
        //p_sample.refreshValidProbes();
        p_sample.setReferenceSegmentForEachProbSet();

        for (Short i : p_Collection.ChrIdInterested) {
            //debug information
            p_Collection.printSegmentsString(curSampleId, i);
        }
        //output Half copy number histogram of segments
        Histogram h = new Histogram(128, 0, 5);
        for (Short i : p_Collection.ChrIdInterested) {
            //debug information
            Set<Segment> segments = p_Collection.getChromosome(curSampleId, i).segments;
            for (Segment seg : segments) {
                h.addDataPoint(seg.robustMeanHalfCopyNumber, seg.Loci.size());
            }
        }
        m_log.info(h.toString("Segments"));

        minCN_segment = 0;
        for (int i = 1; i < h.getCount(); i++) {
            if (h.getFreq(i) > 0) {
                minCN_segment = (h.getThresh(i - 1) + h.getThresh(i)) / 2;
                break;
            }
        }
        m_log.info(String.format("minimum CN of segments = %.4f", minCN_segment));

        CNPeakArray = getPeaks(h);

    }


    /**
     * Clear cutter.
     */
    public void clearCutter() {
        //m_cutter = null;
        for (short i : p_Collection.getChrIdArray()) {
            for (Segment s : p_sample.getAllSegments(i)) {
                s.releaseIntegralCN_Index();
            }
        }

    }

    /**
     * Normalize by mean.
     */
    public void normalizeByMean() {
        double sum = 0;
        int count = 0;
        for (Short index : p_Collection.ChrIdInterested) {
            for (ProbeSet p : p_sample.getAllProbeSets(index).values()) {
                if (p.data.Is_Valid) {
                    sum += p.data.CopyNumberHalf;
                    count++;
                }
            }
        }
        double avg = sum / count;
        for (Short index : p_Collection.ChrIdInterested) {
            for (ProbeSet p : p_sample.getAllProbeSets(index).values()) {
                if (p.data.Is_Valid) {
                    p.data.CopyNumberHalf /= avg;
                }
            }
        }

    }

    private void correctCrossTalk() {
        List<Double> CrosstalkListNormal = new LinkedList<Double>();

        List<Double> CrosstalkListTumor = new LinkedList<Double>();

        for (short i : p_Collection.getChrIdArray()) {
            for (ProbeSet ps : p_sample.getAllProbeSets(i).values()) {
                if (ps.Is_SNPProbe && ps.data.Is_Valid) {
                    ProbeSetRecord nPr = ps.data.Normal;
                    ProbeSetRecord tPr = ps.data.Tumor;

                    if (nPr.isGenotypeA2B2) {
                        if (nPr.mean_A > 0) CrosstalkListNormal.add(nPr.mean_A);
                        if (nPr.mean_B > 0) CrosstalkListNormal.add(nPr.mean_B);
                    }

                    if (tPr.isGenotypeA2B2) {
                        if (tPr.mean_A > 0) CrosstalkListTumor.add(tPr.mean_A);
                        if (tPr.mean_B > 0) CrosstalkListTumor.add(tPr.mean_B);
                    }
                }
            }
        }
        double meanValueNormal = BioToolbox.mean(CrosstalkListNormal);
        double meanValueTumor = BioToolbox.mean(CrosstalkListTumor);

        double np2 = 1.0 / (1.0 - meanValueNormal);
        double np1 = meanValueNormal / (meanValueNormal - 1);
        double tp2 = 1.0 / (1.0 - meanValueTumor);
        double tp1 = meanValueTumor / (meanValueTumor - 1);

        for (short i : p_Collection.getChrIdArray()) {
            for (ProbeSet ps : p_sample.getAllProbeSets(i).values()) {
                if (ps.Is_SNPProbe && ps.data.Is_Valid) {

                    double an = np1 * ps.data.Normal.mean_A + np2 * ps.data.Normal.mean_B;
                    double bn = np2 * ps.data.Normal.mean_A + np1 * ps.data.Normal.mean_B;
                    double at = tp1 * ps.data.Tumor.mean_A + tp2 * ps.data.Tumor.mean_B;
                    double bt = tp2 * ps.data.Tumor.mean_A + tp1 * ps.data.Tumor.mean_B;
                    double abn = an + bn;

                    if (abn == 0) {
                        ps.data.Is_Valid = false;
                    } else {
                        ps.data.alleleA_CorrectCrossTalk = at / abn;
                        ps.data.alleleB_CorrectCrossTalk = bt / abn;
                    }
                }
            }
        }
    }

    private boolean isCalculatedForRou(ProbeSet p) {
        return (p.Is_SNPProbe) && (p.data.Normal.isGenotypeAB) && p.data.Is_Valid;
//        return (p.Is_SNPProbe) && p.data.Is_Valid; //Only one peak
//                    if (((p.Is_SNPProbe) && (!p.data.Normal.isGenotypeA2B2) && p.data.Is_Valid)) {
    }

    /**
     * Calculate correlation coefficient.
     *
     * @param WindowsSize the windows size
     */
    public void calculateCorrelationCoefficient(int WindowsSize) {
        if (p_config.params.enableCorrectCrossTalk) {
            correctCrossTalk();
        }

        int width = WindowsSize | 0x01;
        int halfwidth = width / 2;
        PearsonsCorrelation pCorr = new PearsonsCorrelation();

        for (Short i : p_Collection.ChrIdInterested) {
            for (Segment s : p_sample.getAllSegments(i)) {
                Map<Long, ProbeSet> mps = p_sample.getAllProbeSets(i);

                List<Double> AlleleA = new ArrayList<Double>();
                List<Double> AlleleB = new ArrayList<Double>();

                List<ProbeSet> psList = new LinkedList<ProbeSet>();

                for (Long id : s.Loci.keySet()) {
                    ProbeSet p = mps.get(id);

                    if (isCalculatedForRou(p)) {
                        if (p_config.params.enableCorrectCrossTalk) {
                            AlleleA.add(p.data.alleleA_CorrectCrossTalk);
                            AlleleB.add(p.data.alleleB_CorrectCrossTalk);
                        } else {
                            AlleleA.add(p.data.alleleA);
                            AlleleB.add(p.data.alleleB);
                        }
                        psList.add(p);
                    }
                }

                int len = psList.size();
                if (len < WindowsSize) {
                    continue;
                }


                double[] a = new double[len];
                double[] b = new double[len];
                for (int j = 0; j < len; j++) {
                    a[j] = AlleleA.get(j);
                    b[j] = AlleleB.get(j);
                }

//                a = BioToolbox.blur(a, 10);
//                b = BioToolbox.blur(b, 10);
//                a = BioToolbox.GaussianBlur(a, 10, 1.0);
//                b = BioToolbox.GaussianBlur(b, 10, 1.0);

                double[] aa, bb;

                Iterator<ProbeSet> it = psList.iterator();
                for (int k = 0; k < len; k++) {
                    ProbeSet p = it.next();
                    int start = Math.max(0, k - halfwidth);
                    int stop = Math.min(len, k + halfwidth);

                    aa = Arrays.copyOfRange(a, start, stop);
                    bb = Arrays.copyOfRange(b, start, stop);
                    p.data.rouIsValid = true;
                    p.data.rou = pCorr.correlation(aa, bb);
                }
            }
        }
    }

    /**
     * Find allele balanced probe.
     */
    public boolean findAlleleBalancedProbe() {
        int width = 128;

        Histogram h = new Histogram(width, -1.0, 1.0);

        for (Chromosome chromosome : p_sample.getAllChromosomes()) {
            for (ProbeSet p : chromosome.probes.values()) {
                if (isCalculatedForRou(p) && p.data.rouIsValid) {
                    h.addDataPoint(p.data.rou);
                    // 还需要改进,已经在下面用了新的思路
                    //p.data.isAlleleBalanced = (Math.abs(p.data.rou) < 0.3);
                }
            }
        }

        LOHArray = getPeaks_Valleys(h, 21);

        double[] LOHThresh = new double[2];
        double rightPeak = 1;
        int peakCount = 0;
        LOHThresh[0] = -1;
        LOHThresh[1] = 1;

        for (Map.Entry<Double, Boolean> point : LOHArray.descendingMap().entrySet()) {
            double pv = point.getKey();
            boolean isPeak = point.getValue();
            if (isPeak) {
                if (peakCount == 0) {
                    rightPeak = pv;
                }
                peakCount++;
            } else {
                if (peakCount == 1) { //get first valley from right to left
                    LOHThresh[0] = pv;
                }
            }
        }

        double percent = h.percent(LOHThresh[0]);

        if (peakCount > 0) {
            if ((percent < 0.3)) {
                m_log.severe("\n LOHThresh move to the right peak.");

                if ((rightPeak > -0.3)) {
                    LOHThresh[0] = rightPeak;
//                    peakCount = 2;
                } else {
                    peakCount = 0;
                }
            } else if (Math.abs(LOHThresh[0]) > 0.3) {
                LOHThresh[0] = Math.signum(LOHThresh[0]) * 0.3;
            }
        }

        if (peakCount > 1) {

            for (Chromosome chromosome : p_sample.getAllChromosomes()) {
                for (ProbeSet p : chromosome.probes.values()) {
                    if (p.data.rouIsValid &&
                            (p.data.rou > LOHThresh[0]) && p.data.rou <= LOHThresh[1]) {
                        p.data.isNotLOH = true;
                    }
                }
            }

            m_log.info(h.toString("rou Histogram of AB-Genotype Loci using Half CopyNumber"));

            if (p_config.params.enable_rou_chart) {
                String temp = String.format("LOH Thresh = %.4f", LOHThresh[0]);
                h.saveToPicture("CN_AB_Rou", curSampleId, temp, "rou", "count");
            }
            m_log.info(String.format("LOH Thresh = [%.4f, %.4f]", LOHThresh[0], LOHThresh[1]));

            return true;
        } else {
            m_log.severe("\n don't find the proper LOHThresh");

            if (p_config.params.enable_rou_chart) {
                m_log.info(h.toString("rou Histogram of AB-Genotype Loci using Half CopyNumber"));
                h.saveToPicture("CN_AB_Rou", curSampleId, "LOH is NULL", "rou", "count");
            }

            m_log.info(String.format("LOH Thresh is NULL, Peak count = %d", peakCount));

            return false;
        }

    }

    public void scanCorrelationCoefficient(int WindowsSize) {
        if (p_config.params.enableCorrectCrossTalk) {
            correctCrossTalk();
        }

        PearsonsCorrelation pCorr = new PearsonsCorrelation();

        Histogram[] hArray = new Histogram[6];
        for (int n = 0; n < 6; n++) {
            hArray[n] = new Histogram(128, -1.0, 1.0);
        }

        for (Short i : p_Collection.ChrIdInterested) {
            for (Segment s : p_sample.getAllSegments(i)) {
                Map<Long, ProbeSet> mps = p_sample.getAllProbeSets(i);

                List<Double> AlleleA = new ArrayList<Double>();
                List<Double> AlleleB = new ArrayList<Double>();

                for (Long id : s.Loci.keySet()) {
                    ProbeSet p = mps.get(id);

                    if (isCalculatedForRou(p)) {
                        if (p_config.params.enableCorrectCrossTalk) {
                            AlleleA.add(p.data.alleleA_CorrectCrossTalk);
                            AlleleB.add(p.data.alleleB_CorrectCrossTalk);
                        } else {
                            AlleleA.add(p.data.alleleA);
                            AlleleB.add(p.data.alleleB);
                        }
                    }
                }

                int len = AlleleA.size();

                double[] a = new double[len];
                double[] b = new double[len];
                for (int j = 0; j < len; j++) {
                    a[j] = AlleleA.get(j);
                    b[j] = AlleleB.get(j);
                }

                double[] aa, bb;

                for (int j = 0; j < 5; j++) {
                    int width = ((WindowsSize >> 2) << j) | 0x01;
                    int halfwidth = width / 2;

                    if (len < width) {
                        continue;
                    }

                    for (int k = 0; k < len; k++) {
                        int start = Math.max(0, k - halfwidth);
                        int stop = Math.min(len, k + halfwidth);

                        aa = Arrays.copyOfRange(a, start, stop);
                        bb = Arrays.copyOfRange(b, start, stop);
                        double c = pCorr.correlation(aa, bb);
                        hArray[j].addDataPoint(c);
                        hArray[5].addDataPoint(c);
                    }
                }
            }
        }

        ChartXYLine cxy = new ChartXYLine();

        for (int j = 0; j < 6; j++) {
            XYSeries xyseries = new XYSeries("Rou" + Integer.toString(j));

            if (j < 5) {
                for (int k = 1; k < hArray[j].getBin(); k++) {
                    xyseries.add(hArray[j].getThresh(k), hArray[j].getFreq(k));
                }
            } else {
                for (int k = 1; k < hArray[j].getBin(); k++) {
                    xyseries.add(hArray[5].getThresh(k), hArray[5].getFreq(k) / 5.0);
                }

            }

            cxy.addLine(xyseries);
        }
        cxy.createChart("Rou", "x", "y");
        XYPlot xyPlot = (XYPlot) cxy.jfreechart.getPlot();
        XYLineAndShapeRenderer lineandshaperenderer = (XYLineAndShapeRenderer) xyPlot.getRenderer();
        lineandshaperenderer.setSeriesStroke(5, new BasicStroke(3.0F));

        Date date = new Date();
        SimpleDateFormat df = new SimpleDateFormat("yyyyMMdd_HHmmss");
        String filename = String.format("..%1$sData%1$s%2$s_%3$s_%4$s.png",
                File.separator, "ScanRou", p_sample.getSampleId(), df.format(date));

        cxy.saveAsFile(filename, 600, 400);

    }

    private TreeMap<Integer, Double> getPeaks(Histogram hist) {
        TreeMap<Integer, Double> peakMap = new TreeMap<Integer, Double>();

        int[] Int_freq = hist.getFreqArray();
        double[] freq = BioToolbox.GaussianBlur(Int_freq, 5, 1);
        double mean = BioToolbox.meanEx(freq);

        //求差分
        double[] diffcn = new double[freq.length];
        diffcn[0] = 0;
        for (int j = 1; j < freq.length; j++) {
            if (freq[j] > mean) {
                diffcn[j] = (freq[j] - freq[j - 1]);
            }
        }
        diffcn = BioToolbox.GaussianBlur(diffcn, 3, 1);

        for (int j = 2; j < diffcn.length - 2; j++) {
            if ((diffcn[j] > 0) && (diffcn[j + 1] <= 0)) {
                peakMap.put(new Integer((int) freq[j]), (hist.getThresh(j - 1) + hist.getThresh(j)) / 2);
            }
        }
        return peakMap;
    }

    private TreeMap<Double, Boolean> getPeaks_Valleys(Histogram hist, int width) {
        TreeMap<Double, Boolean> peakMap = new TreeMap<Double, Boolean>();

        int[] Int_freq = hist.getFreqArray();
        double[] freq = BioToolbox.GaussianBlur(Int_freq, width, 1);
        double mean = BioToolbox.meanEx(freq);

        //求差分
        double[] diffcn = new double[freq.length];
        diffcn[0] = 0;
        for (int j = 1; j < freq.length; j++) {
//            if (freq[j] > mean) {
            diffcn[j] = (freq[j] - freq[j - 1]);
//            }
        }
        diffcn = BioToolbox.GaussianBlur(diffcn, 3, 1);

        for (int j = 2; j < diffcn.length - 2; j++) {
            if ((diffcn[j] > 0) && (diffcn[j + 1] <= 0)) {
                peakMap.put((hist.getThresh(j - 1) + hist.getThresh(j)) / 2, true); //is peak
            } else if ((diffcn[j] < 0) && (diffcn[j + 1] >= 0)) {
                peakMap.put((hist.getThresh(j - 1) + hist.getThresh(j)) / 2, false); //is valley
            }
        }
        return peakMap;
    }

    /**
     * Find allele balance peak.
     *
     * @return the double [ ]
     */
    public double[] findAlleleBalancePeak() {
        final int count = 128;
        Histogram h = new Histogram(count, 0, 4);

        for (Short i : p_Collection.getChrIdArray()) {
            Map<Long, ProbeSet> chromosome = p_sample.getAllProbeSets(i);
            for (Segment s : p_sample.getAllSegments(i)) {
                double avg = s.robustMeanHalfCopyNumber;
                for (Long id : s.Loci.keySet()) {
                    ProbeSet ps = chromosome.get(id);
                    if (isCalculatedForRou(ps) && ps.data.isNotLOH) {
                        h.addDataPoint(avg);
                    }
                }
            }
        }
        m_log.info(h.toString("AB Normal part of copynumber"));

        TreeMap<Integer, Double> peakMap = getPeaks(h);

        double peakValue = peakMap.lastEntry().getValue();

        //There are two peaks of balanced Loci.
        if (peakMap.size() > 1) {
            double secondValue = peakMap.lowerEntry(peakMap.lastKey()).getValue();

            double fourthPeak = Math.max(peakValue, secondValue);
            peakValue = Math.min(peakValue, secondValue);

            if (fourthPeak > 1.15 * peakValue) { //There are two peaks
                p_sample.fourfoldPeak = fourthPeak;
                m_log.info(String.format("%n\t\t\t\tEstimated Peak of (CN=4.0) = %.4f", fourthPeak));

                double r4 = testingCopyNumberEqualFour(fourthPeak);
                double r2 = testingCopyNumberEqualFour(peakValue);
//                System.exit(0);

                p_sample.getChromosome(0).baseline = peakValue;
                return new double[]{peakValue * (1 - p_config.params.widthNormalPart),
                        peakValue * (1 + p_config.params.widthNormalPart)};
            } else {//too close
                peakValue = peakMap.lastEntry().getValue();
            }
        }

        //maybe peakValue refer to copynumber = 4.0 or 2.0
        int n = 0;
        CNPeakArray.keySet();
        for (double peak : CNPeakArray.values()) {
            if (peakValue - peak > 0.1) {
                n++;
                if (n == 2) {
                    p_sample.fourfoldPeak = peakValue;
                    peakValue = peak;
                    p_sample.getChromosome(0).baseline = peakValue;
                    return new double[]{peakValue * (1 - p_config.params.widthNormalPart),
                            peakValue * (1 + p_config.params.widthNormalPart)};
                }
            }
        }

        double r4 = testingCopyNumberEqualFour(peakValue);
        m_log.info(String.format("# factor of 4N = %.3f", r4));
        if (r4 < 3) { // copy number == 4
            p_sample.fourfoldPeak = peakValue;
            peakValue /= 2;
        } else if (r4 > 7) { // copy nubmer == 2

        } else { // maybe 4 or 2

        }
        p_sample.getChromosome(0).baseline = peakValue;
        return new double[]{peakValue * (1 - p_config.params.widthNormalPart),
                peakValue * (1 + p_config.params.widthNormalPart)};

    }

    private double testingCopyNumberEqualFour(double peakValue) {
        double down = peakValue * (1 - p_config.params.widthNormalPart);
        double up = peakValue * (1 + p_config.params.widthNormalPart);
        Histogram h = new Histogram(64, 0, 1);
        for (Short i : p_Collection.getChrIdArray()) {
            Map<Long, ProbeSet> chromosome = p_sample.getAllProbeSets(i);
            for (Segment s : p_sample.getAllSegments(i)) {
                double avg = s.robustMeanHalfCopyNumber;
                if (avg < up && avg > down) {
                    for (Long id : s.Loci.keySet()) {
                        ProbeSet ps = chromosome.get(id);
                        if (ps.data.Is_Valid && ps.Is_SNPProbe) {
                            h.addDataPoint(ps.data.BAF);
                        }
                    }
                }
            }
        }

        double[] dBAF = new double[h.getBin() + 1];
        int sum = h.getCount();
        for (int i = 0; i < h.getBin(); i++) {
            dBAF[i] = ((double) h.getFreq(i + 1)) / sum;
        }
        FastCosineTransformer dct = new FastCosineTransformer(DctNormalization.ORTHOGONAL_DCT_I);
        double[] dctBAF = Arrays.copyOfRange(dct.transform(dBAF, TransformType.FORWARD), 1, dBAF.length);
        //dctBAF = BioToolbox.GaussianBlur(dctBAF,3,1);
        double tag10 = 0, tag20 = 0;
        for (int i = 0; i < 20; i++) {
            if (i < 10) {
                tag10 += Math.abs(dctBAF[i]);
            } else {
                tag20 += Math.abs(dctBAF[i]);
            }
        }
//        for (int i = 0; i < dBAF.length; i++){
//            System.out.println(dBAF[i]);
//        }
//        System.out.println("##############################");

        return (tag10 / tag20);
    }


    /**
     * Check deletion.
     */
    public void checkDeletion() {

        double[] range = findAlleleBalancePeak();

        m_log.info(String.format("Range of Half copynumber of copy-neutral Loci: [%.5f - %.5f]", range[0], range[1]));

        for (Short i : p_Collection.getChrIdArray()) {
            //Map<Long, ProbeSet> chromosome = p_Collection.sample.get(i).probes;
            for (Segment s : p_sample.getAllSegments(i)) {
                if (s.fVector == null) {
                    s.fVector = new SegmentFeature();
                    s.fVector.Seg_id = s.Seg_id;
                }

                if ((s.robustMeanHalfCopyNumber >= range[0]) && (s.robustMeanHalfCopyNumber <= range[1])) {
                    s.fVector.isNeutral = true;
                    s.fVector.isDeletion = false;
                } else if (s.robustMeanHalfCopyNumber < range[0]) {
                    s.fVector.isDeletion = true;
                    s.fVector.isNeutral = false;
                } else {
                    s.fVector.isNeutral = false;
                    s.fVector.isDeletion = false;
                }
            }
        }
    }

    /**
     * Calcuate base line.
     */
    public void calcuateBaseLine() {
        double sumChr, sumAll;
        int countChr, countAll;

        //Histogram h = new Histogram(128, 0, 8);
        sumAll = 0;
        countAll = 0;
        for (Short i : p_Collection.getChrIdArray()) {
            Map<Long, ProbeSet> chromosome = p_sample.getAllProbeSets(i);
            sumChr = 0;
            countChr = 0;

            for (Segment s : p_sample.getAllSegments(i)) {
                double segCN = s.robustMeanHalfCopyNumber;
                if (s.fVector.isNeutral) {
                    for (Long id : s.Loci.keySet()) {
                        ProbeSet ps = chromosome.get(id);
                        //条件值得考虑，是否加入Normal.isGenotypeAB好像没有什么影响
                        if ((ps.data.Is_Valid) && (ps.data.Normal.isGenotypeAB)) {
                            sumChr += segCN;
                            countChr++;
                        }
                    }
                }
            }

            if (countChr > p_config.params.MinimumCountForBaseLine) {
                p_sample.getChromosome(i).baseline = sumChr / countChr;
            } else {
                p_sample.getChromosome(i).baseline = null;
            }
            sumAll += sumChr;
            countAll += countChr;
        }
        if (countAll > 0) {
            p_sample.getChromosome(0).baseline = sumAll / countAll;
        }

    }

    /**
     * Normalize by base line.
     */
    public void normalizeByBaseLine() {
        double baseline;
        //baseline = p_sample.getChromosome(0).baseline;

        for (int i : p_Collection.getChrIdArray()) {
            Map<Long, ProbeSet> chromosome = p_sample.getAllProbeSets(i);

            //baseline = p_sample.getChromosome(0).baseline;
            if (p_sample.getChromosome(i).baseline != null) {
                baseline = p_sample.getChromosome(i).baseline;
            } else {
                baseline = p_sample.getChromosome(0).baseline;
            }

            for (ProbeSet p : chromosome.values()) {
                if (p.data.Is_Valid) {
                    p.data.CopyNumberHalf = p.data.CopyNumberHalf / baseline;
                    if (p.Is_SNPProbe) {
                        p.data.alleleA = p.data.alleleA / baseline;
                        p.data.alleleB = p.data.alleleB / baseline;
                        p.data.LRR = BioToolbox.log2(2 * (p.data.alleleA + p.data.alleleB));
                    }
                }
            }

            for (Segment s : p_sample.getAllSegments(i)) {
                s.avgHalfCopyNumber = s.avgHalfCopyNumber / baseline;
                s.maxHalfCopyNumber = s.maxHalfCopyNumber / baseline;
                s.minHalfCopyNumber = s.minHalfCopyNumber / baseline;
                s.robustMeanHalfCopyNumber = s.robustMeanHalfCopyNumber / baseline;
            }
        }

        p_sample.fourfoldPeak /= p_sample.getChromosome(0).baseline;
    }

    /**
     * Calc deletion histogram.
     */
    public void calcDeletionHistogram() {
        Histogram h = new Histogram(32, 0, 2);
        for (int i : p_Collection.getChrIdArray()) {
            // Map<Long, ProbeSet> chromosome = p_sample.getAllProbeSets(i);
            for (Segment s : p_sample.getAllSegments(i)) {
                if (s.fVector.isDeletion) {
//                    for (Long id : s.Loci.keySet()){
//                        ProbeSet ps = chromosome.get(id);
//                        if (ps.data.Is_Valid){
//                            h.addDataPoint(ps.data.CopyNumberHalf);
//                        }
//                    }
                    h.addDataPoint(s.robustMeanHalfCopyNumber, s.Loci.size());
                }

            }
        }
        m_log.info(h.toString("Deletion Segment"));
    }
}
//double[] range = findAlleleBalancePeak();
//
//        double chrSum = 0;
//        double sampleSum = 0;
//        int chrCount = 0;
//        int sampleCount = 0;
//
//        for (int i = 1; i < 27; i++) {
//            chrSum = 0;
//            chrCount = 0;
//            Map<Long, ProbeSet> chromosome = p_Collection.sample.get(i).probes;
//
//            for (Segment s : p_Collection.sample.get(i).segments) {
//                if (s.fVector == null) {
//                    s.fVector = new SegmentFeature();
//                    s.fVector.Seg_id = s.Seg_id;
//                }
//
//                if ((s.robustMeanHalfCopyNumber > range[0]) && (s.robustMeanHalfCopyNumber < range[1])) {
//                    s.fVector.isNeutral = true;
//                    s.fVector.isDeletion = false;
//
//                    for (Long id : s.Loci.keySet()) {
//                        ProbeSet ps = chromosome.get(id);
//                        ps.data.fVector.Seg_id = s.Seg_id;
//                        //if ((ps.data.Is_Valid)&&(ps.data.fVector.isAlleleBalanced) && (ps.data.fVector.isGenotypeAB)) {
//                        if (ps.data.Is_Valid) {
//                            chrSum += ps.data.CopyNumberHalf.copynumber;
//                            chrCount++;
//                        }
//                    }
//                } else if (s.avgHalfCopyNumber <= range[0]) {
//                    s.fVector.isDeletion = true;
//                    s.fVector.isNeutral = false;
//                    //System.out.println("Deletion:"+ s.avgHalfCopyNumber);
//                } else {
//                    s.fVector.isNeutral = false;
//                    s.fVector.isDeletion = false;
//                }
//            }
//
//            if (chrCount > 0) {
//                p_Collection.sample.get(i).baseline = chrSum / chrCount;
//
//                System.out.println(String.format("chr_id = %2d;\t count = %d\t; baseline = %f \t [%.5f - %.5f]", i, chrCount,
//                        p_Collection.sample.get(i).baseline, range[0], range[1]));
//
//            } else {
//                p_Collection.sample.get(i).baseline = null;
//            }
//
//            sampleSum += chrSum;
//            sampleCount += chrCount;
//        }
//        p_Collection.sample.get(0).baseline = sampleSum / sampleCount;


//        List<Double> peaks = new ArrayList<Double>();
//        for (int j = 2; j < diffcn.length - 2; j++) {
//            //if (((diffcn[j]<0)&&(diffcn[j+1]>0)) || ((diffcn[j]>0)&&(diffcn[j+1]<0))){
//            if ((diffcn[j] < 0) && (diffcn[j + 1] > 0)) {//波谷
//                peaks.add(h.getThresh(j - 1));
//            }
//        }


//        double LOHThresh;
//        if (peaks.size() > 0){
//            LOHThresh = peaks.get(peaks.size()-1);
//        }else{
//            LOHThresh = -0.2;
//        };

//    public void combineSegments() {
//
//        for (Short i : p_Collection.getChrIdArray()) {
//            Set<Segment> newSegs = new TreeSet<Segment>();
//            Segment priorSeg;
//            priorSeg = null;
//            for (Segment s : p_Collection.sample.get(i).segments) {
//                if (newSegs.size() > 0) {
//                    if ((priorSeg != null) && s.compareSegmentKind(priorSeg)) {
//                        //combine
//                        priorSeg.Loci.putAll(s.Loci);
//                    } else {
//                        newSegs.add(s);
//                        priorSeg = s;
//                    }
//                } else {
//                    newSegs.add(s);
//                    priorSeg = s;
//                }
//            }
//
//            for (Segment ns : newSegs) {
//                ns.updateCharacters();
//            }
//            p_Collection.sample.get(i).segments = newSegs;
//        }
//    }


//    public void calcuateBaseLine() {//计算局部BaseLine
//        List<Double> oneChr = new ArrayList<Double>();
//        List<Double> base = new ArrayList<Double>();
//
//        for (Short i : p_Collection.getChrIdArray()) {
//            Map<Long, ProbeSet> chromosome = p_sample.getAllProbeSets(i);
//            oneChr.clear();
//
//            for (Segment s : p_sample.getAllSegments(i)) {
//                if (s.fVector.isNeutral) {
//                    for (Long id : s.Loci.keySet()) {
//                        ProbeSet ps = chromosome.get(id);
//                        //条件值得考虑，是否加入Normal.isGenotypeAB好像没有什么影响
//                        if ((ps.data.Is_Valid) && (ps.data.Normal.isGenotypeAB)) {
//                            oneChr.add(ps.data.CopyNumberHalf);
//                        }
//                    }
//                }
//            }
//            if (oneChr.size() > 500) {
//                double[] values = BioToolbox.robustMean(oneChr, 16);
//
//                p_sample.getChromosome(i).baseline = values[0];
//                p_sample.getChromosome(i).baseline_std = values[1];
//                base.add(values[0]);
//            } else {
//                p_sample.getChromosome(i).baseline = null;
//            }
//        }
//
//        p_sample.getChromosome(0).baseline = BioToolbox.mean(base);
//    }
//    public void checkDeletion(){
//        for (Short i : p_Collection.getChrIdArray()){
//            Map<Long, ProbeSet> chromosome = p_Collection.sample.get(i).probes;
//            for (Segment s : p_Collection.sample.get(i).segments) {
//                if (s.fVector == null) {
//                    s.fVector = new SegmentFeature();
//                    s.fVector.Seg_id = s.Seg_id;
//                }
//
//                if ((s.robustMeanHalfCopyNumber > 0.8) && (s.robustMeanHalfCopyNumber < 1.2)) {
//                    s.fVector.isNeutral = true;
//                    s.fVector.isDeletion = false;
//                } else if (s.avgHalfCopyNumber <= 0.8) {
//                    s.fVector.isDeletion = true;
//                    s.fVector.isNeutral = false;
//                } else {
//                    s.fVector.isNeutral = false;
//                    s.fVector.isDeletion = false;
//                }
//            }
//        }
//    }
//方案三：零点距离，负的
//            int begin = valleys.size() - 1;
//            LOHThresh[0] = h.getThresh(valleys.get(begin));
//            for (int i = begin - 1; i > 0; i--) {
//                double tt = h.getThresh(valleys.get(i));
//                if (Math.abs(tt) < Math.abs(LOHThresh[0])) {
//                    LOHThresh[0] = tt;
//                }
//            }

//    public double[] findAlleleBalancePeak2() {
//        Histogram h = new Histogram(128, 0, 5);
//
//        for (Short i : p_Collection.getChrIdArray()) {
//            Map<Long, ProbeSet> chromosome = p_sample.getAllProbeSets(i);
//            for (Segment s : p_sample.getAllSegments(i)) {
//                double avg = s.robustMeanHalfCopyNumber;
//                for (Long id : s.Loci.keySet()) {
//                    ProbeSet ps = chromosome.get(id);
//                    if (ps.Is_SNPProbe && ps.data.isNotLOH &&
//                            (ps.data.Normal.isGenotypeAB)) {
////                        (ps.data.Normal.isGenotypeAB&& !ps.data.Tumor.isGenotypeA2B2)) {
//                        // or using the mean of copynumber of each segment
//                        h.addDataPoint(avg);
//                    }
//                }
//            }
//        }
//
//        m_log.info(h.toString("AB Normal part of copynumber 1"));
//
//        int[] Int_freq = h.getFreqArray();
//        double[] freq = BioToolbox.GaussianBlur(Int_freq, 5, 1);
//
//        //求差分
//        double[] diffcn = new double[freq.length];
//        diffcn[0] = 0;
//        for (int j = 1; j < freq.length; j++) {
//            diffcn[j] = (freq[j] - freq[j - 1]);
//        }
//        diffcn = BioToolbox.GaussianBlur(diffcn, 3, 1);
//
//        int maxIndex = 0;
//        int secondIndex = 0;
//        for (int j = 2; j < diffcn.length - 2; j++) {
//            if ((diffcn[j] > 0) && (diffcn[j + 1] <= 0)) {
//                if (freq[maxIndex]<freq[j]) {
//                    secondIndex = maxIndex;
//                    maxIndex = j;
//                }else if (freq[secondIndex]<freq[j]){
//                    secondIndex = j;
//                }
//            }
//        }
//
//        double peakValue = (h.getThresh(maxIndex - 1) + h.getThresh(maxIndex)) / 2;
//        double secondValue;
//        if (secondIndex > 0) {
//            secondValue = (h.getThresh(secondIndex - 1) + h.getThresh(secondIndex)) / 2;
//            double secMaxDistance = peakValue - secondValue;
//            double secMinDistance = secondValue - minCN_segment;
//
//            calculatePeakPosition(minCN_segment, peakValue, secondValue);
//
//            //如何引入minCN_segment，增强约束？
//            if ((peakValue > secondValue) && (secMinDistance > secMaxDistance/3)) {
//                return new double[]{secondValue * 0.85, secondValue * 1.15};
//            }
//        }
//
//        return new double[]{peakValue * 0.85, peakValue * 1.15};
//
//    }
//
//    private void calculatePeakPosition(double minValue, double maxPeak, double secondPeak){
//        double step = Math.abs(maxPeak - secondPeak)/3.0;
//        double minStep = Math.abs(Math.min(maxPeak,secondPeak)-minValue);
//
//        if (maxPeak > secondPeak){
//            p_sample.fourthPeakRatio = maxPeak/secondPeak;
//        }else {
//            p_sample.fourthPeakRatio = secondPeak/maxPeak;
//        }
//
//    }